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We show that the topological Kitaev spin liquid on the honeycomb lattice is extremely fragile against the 
second neighbor Kitaev coupling K 2 , which has been recently shown to be the dominant perturbation away 
from the nearest neighbor model in iridate Na 2 lrC> 3 , and may also play a role in a-RuCR and LFIrCH This 
coupling explains naturally the zig-zag ordering (without introducing unrealistically large longer-range Heisen¬ 
berg exchange terms), and the special entanglement between real and spin space observed recently in NaalrOs. 
Moreover, the minimal K 1 -K 2 model that we present here holds the unique property that the classical and 
quantum phase diagrams and their respective order by disorder mechanisms are qualitatively different due to the 
fundamentally different symmetries of the classical and quantum counterparts. 


I. Introduction 

The search for novel quantum states of matter arising from 
the interplay of strong electronic correlations, spin-orbit cou¬ 
pling (SOC), and crystal field splitting has recently gained 
strong impetus in the context of 4 d and 5 d transition metal 
oxides [1], The layered iridates of the A 2 Ir 03 (A=Na,Li) fam¬ 
ily [2-7] have been at the center of this search because of the 
prediction [8, 9] that the dominant interactions in these mag¬ 
nets constitute the celebrated Kitaev model on the honeycomb 
lattice, one of the few exactly solvable models hosting gapped 
and gapless quantum spin liquids (QSLs) [10]. This aspect to¬ 
gether with the realization that the Kitaev spin liquid is stable 
with respect to moderate Heisenberg-like perturbations [9, 11] 
has triggered a lot of experimental activity on A 2 Ir 03 and, 
more recently, on the similar o-RuCI.- compound [12-14], 

In the layered A 2 Ir 03 magnets, the single-ion ground state 
configuration of Ir 4+ is an effective pseudospin J e g = 1/2 
doublet, where spin and orbital angular momenta are in¬ 
tertwined due to the strong SOC. In the original Kitaev- 
Heisenberg model proposed by Jackeli and Khaliullin [8], the 
pseudospins couple via two competing nearest neighbor (NN) 
interactions: An isotropic antiferromagnetic (AFM) Heisen¬ 
berg exchange, Ji, and a highly anisotropic Kitaev interac¬ 
tion, K 1 , which is strong and ferromagnetic, a fact that is 
also confirmed by ab-initio quantum chemistry calculations 
by Katukuri et al [15, 16]. Nevertheless, neither Na 2 Ir 03 nor 
Li 2 Ir 03 are found to be in the spin liquid state at low temper¬ 
atures. Instead, they show, respectively, AFM zigzag and in¬ 
commensurate long-range magnetic orders, none of which is 
actually present in the Kitaev-Heisenberg model for FM K\ 
coupling. 

The most natural way to obtain these magnetic states is 
by including further neighbor Heisenberg couplings [15-18], 
which are non-negligible due to extended nature of the 5 d- 
orbitals of Ir 4+ ions [6, 19]. In addition, recent calculations 
by Sizyuk et al [20] based on the ab-initio density-functional 
data of Foyevtsova et al [21] have shown that, forNa 2 IrC> 3 , the 
next nearest neighbor (NNN) exchange paths must also give 


rise to an anisotropic, Kitaev-like coupling /\ 2 , which turns 
out to be AFM. More importantly, this coupling is the largest 
interaction after I\ 1 . It has also been argued [22] that /\ 2 plays 
an important role in the stabilization of the IC spiral state in 
Li 2 IrC >3 and might be deduced from the strong-coupling limit 
of Hubbard model with topological band structure [23, 24]. 

Recent structural [12] and magnetic [13] studies have 
shown that the layered honeycomb magnet a-RuCl 3 is another 
example of a strong SOC Mott insulator, where the Ru 3+ ions 
are again described by effective J e g = 1/2 doublets. At low 
T, this magnet exhibits zigzag ordering as in Na 2 Ir 03 . Fur¬ 
thermore, the superexchange derivations [25, 26] based on the 
ab initio tight-binding parameters show that the NNN cou¬ 
pling K'z is again appreciable, and the signs of both K\ and 
K -2 are reversed compared to Na 2 Ir 03 (i.e., K\ is AFM and 
K 2 is FM). However, a strong off-diagonal symmetric NN ex¬ 
change F term [15, 16, 27], which is allowed by symmetry, 
is also present [25, 26], together with a much smaller Ji cou- 



FIG. 1. (Color online) The Kitaev K 1 -K 2 model with three types 
of NN (solid) and NNN (dashed) Ising bonds. Here ti = ay and 
t 2 = (—^x+ |y)a are two primitive translations and a is a lat¬ 
tice constant. We also show the vertical 2-leg ladders (shaded strips) 
discussed in the text, and the four-sublattice decomposition (A-D) 
related to the operations H yzx and H xyz , see text. 
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pling. This compound must then be examined in connection 
to T, K 2 , and Ji, since the T term alone is not sufficient to 
explain the experimental situation, as we discuss at length in 
Sec. VII. 

Motivated by these studies, here we consider the minimal 
extension of the NN Kitaev model that incorporates the ef¬ 
fect of K 2 , the K\ -1\‘2 model. We show that an extremely 
weak K 2 is enough to stabilize the zig-zag phases relevant for 
Na 2 Ir 03 and a-RuCI :i , without introducing large, second and 
third neighbor Heisenberg exchange J 2 and J 3 . While J 2 and 
J 3 are present in these compounds, the key point is that the Ki¬ 
taev spin liquid is significantly more fragile against K 2 than 
J 2 and J 3 . Thus, in conjunction with the above predictions 
from superexchange derivations, our findings suggest that any 
adequate minimal model of these compounds should include 
the NNN coupling K 2 . 

A very striking aspect of the zig-zag phases (shared by all 
magnetic phases) of the K\-K 2 model is that they are only 
stabilized for quantum spins and not for classical spins, de¬ 
spite having a strong classical character. Indeed, these phases 
are Ising-like (with spins pointing along one of the three cubic 
axes), they are protected by a large excitation gap in the inter¬ 
acting 1 /S spin-wave spectrum, and the spin lengths are ex¬ 
tremely close to their classical value of 1/2. Yet, these phases 
cannot be stabilized in the classical limit, in stark contrast to 
the conventional situation where quantum and thermal fluctu¬ 
ations work in parallel and often lead to the same order-by- 
disorder phenomena. Instead, this rare situation we encounter 
here stems from the manifestly different symmetry structure 
of the classical and quantum Hamiltonians, and the underlying 
principle that time reversal can only act globally in quantum 
systems (see below). This aspect has important ramifications 
for the phase diagram at zero and finite temperatures T. 

II. Model & Phase diagram 

The model we consider here is described by the effective 
spin-1/2 Hamiltonian 

n = h\J2 s? 1 s] ij + K 2 Y. s i ij S J - (D 

(ij) <ij» 

where (ij) (respectively <^ij)$>) label NN (NNN) spins on the 
honeycomb lattice, 5“ defines the ath cartesian component of 
the spin operator at site j, and 7 ^ (A ij) define the type of Ising 
coupling for the bond (ij), see Fig. 1. This model interpolates 
between two well known limits, the exactly solvable Kitaev 
spin liquid [10] at K 2 = 0, and the triangular Kitaev model 
at A'i = 0 [28-32], It is easy to see that a finite K 2 ruins the 
exact solvability of the NN Kitaev model because the flux op¬ 
erators [10] W p = 2 6 S{S 2 SfS'fS'f Sq (see site-labeling con¬ 
vention in Fig. 5, top left), around hexagons p are no longer 
conserved. 

In the following we parametrize Ki = cos ip and K 2 = 
sin ip, and take ip £ [0, 2n). It turns out that the physics ac¬ 
tually remains the same under a simultaneous sign change of 
Ki and K 2 , because this can be gauged away by an operation 



FIG. 2. (Color online) (a) The T = 0 phase diagram of the model 
(1) as found by exact diagonalizations. Each of the magnetic regions 
(I-IV) hosts twelve degenerate quantum states. Here we show two 
members (where spins point along the z-axis, blue/red circles denote 
spin up/down) that are related to each other by flipping the spins in 
every second ladder (shaded strips) of Fig. 1. The Bragg peaks corre¬ 
sponding to ( SfSj) correlations are also shown in the extended Bril- 
louin zone (assuming the same magnetic form factor in the two unit 
cell sublattices). The corresponding Bragg reflections for (SfSj) 
and (S'/ Sj) are related to (Sf SJ) by Cg v spin-orbit rotations [33], 


Hyzx = n j6 B c 2 y(i) n,ec c 2z(j) rifceD c 2 x(k), which is the 
product of 7r-rotations around the y, z, and x axis, respec¬ 
tively, for the B, C, and D sublattices of Fig. 1. This hidden 
duality is a very common feature in many spin-orbital mod¬ 
els [9, 34, 35] but does not exist when Heisenberg couplings 
are also present (in contrast to the symmetry H xyz discussed 
below). Here it reduces our study to the first two quadrants of 
the unit circle of ip. 

Figure 2 shows the quantum phase diagram as found by 
exact diagonalizations (ED) on finite clusters, see discussion 
below and numerical data shown in Fig. 3. There are six dif¬ 
ferent regimes as a function of the angle ip: the two quan¬ 
tum spin liquids (QSLs) regions (which have been enlarged 
for better visibility) around the exactly solvable Kitaev points 
(ip = 0 and 7 r) and four long-range magnetic regions (I-IV), 
hosting FM, Neel, stripy, as well as the zig-zag phases that 
are relevant for Na 2 IrC >3 (II) and 0-R11CI3 (IV). Under the du¬ 
ality transformation H yzx , the two QSLs map to each other, I 
maps to III, and II maps to IV. 

Each of the magnetic regions actually hosts twelve degen¬ 
erate quantum states, some of which are even qualitatively 
different among themselves, with very distinct Bragg reflec¬ 
tions. For example, the region III hosts six FM and six stripy 
AFM ground states, and IV hosts six Neel and six zigzag AFM 
ground states. This striking aspect stems from a non-global 
symmetry, H xyz , which is the product of 7r-rotations around 
the x, y, and z axis, respectively, for the B, C, and D sublat¬ 
tices of Fig. 1. The two states shown in each magnetic region 
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FIG. 3. (Color online) (a-b) Exact low-energy spectra (measured from the ground state energy Eo) of the 24-site (a) and 32-site (b) clusters, 
defined, respectively [33], by the spanning vectors (2ti—4t2,4ti— 2t2) and (2ti—4t2,4ti). A non-linear x-axis is used in order to highlight 
all regions of interest equally. The states are labeled by momenta k in the first BZ, parity (“e” for even, “o” for odd) under inversion through 
hexagon centers, and parity under global spin 7r-rotations around the x-axis (“Sze” for even, “Szo” for odd). The (red) numbers in (a) denote 
the multiplicity of the lowest five levels in regions I and II, and the ground state degeneracy at ip = 0 and 7r. (c) Ground state expectation value 
( W p ) of Kitaev's flux operators, (d) Square root of the ‘symmetrized" ground state spin structure factor <S(Q) (see text), along with the spin 
length calculated from a self-consistent non-linear spin-wave theory (NLSWT). 


of Fig. 2 are related to each other by this symmetry, which for 
these particular states amounts to flipping the z-component of 
the spins in every second shaded ladder of Fig. 1. The remain¬ 
ing ten states of the quantum ground state manifold arise by 
applying the global symmetries of the model: i) the double 
cover Cg v of C ( , v , and ii) the double cover D 2 of the D 2 group 
of global 7t rotations in spin space. 

Let us now turn to the numerical spectra shown in 
Fig. 3 (a,b). First, the QSL regions are extremely narrow: 
They survive in a tiny window of 5 ip = 0.057r around the exact 
Kitaev points, which is confirmed by the comparison of ED 
against large scale pseudofermion functional renormalization 
group (PFFRG) calculations [36-39], So the QSLs are ex¬ 
tremely fragile against K 2 . 

Second, Fig. 3 (a,b) show very dense spectral features in the 
QSL regions, reflecting the continuum structure of fractional¬ 
ized excitations above the Kitaev spin liquid. More specifi¬ 
cally, for finite systems the ground state degeneracy at the ex¬ 
act Kitaev points [40] is lifted by K 2 . Still, for small enough 
\K 2 \, the QSLs must be gapless in the thermodynamic limit, 
because K 2 respects time reversal symmetry and is therefore 
not expected [10] to open a gap in the Majorana spectrum [41], 

Third, unlike the QSL regions, the low-energy spectrum 
inside the magnetic regions is very discrete. In addition, 
most of the low-lying states within the energy window shown 
in Figs. 3 (a,b) correspond precisely to the twelve quantum 
ground states discussed above. For finite systems, these states 
are admixed by a finite tunneling, leading to twelve symmet¬ 
ric eigenstates with quantum numbers corresponding to the 
decomposition of the symmetry broken states. This decompo¬ 
sition is worked out in detail in [33] and is indeed fully con¬ 
sistent with the ED data. So the lowest twelve states in each 
magnetic region of Figs. 3 (a,b) will collapse to zero energy in 


the thermodynamic limit, leaving the true magnon excitations 
with a large anisotropy gap (modulo finite size corrections), 
reflecting the anisotropic, Ising-like character of the magnetic 
model. 

Fourth, the magnetic instabilities, which serve as good ex¬ 
amples of deconfinement-confinement transitions [42-45] for 
the underlying spinons, are of first order, as they are accom¬ 
panied by finite, abrupt changes [46] in several ground state 
properties, e.g., in ( W p ), and in the spin-spin correlations. 
Specifically, at ip = 0 and n, all fluxes W p have a value of 
+1 [10]. A finite I\ 2 admixes sectors of different W p , and 
so (W p ) drops continuously as we depart from the exact Ki- 
taev’s points, until it jumps to very low absolute values when 
we enter the magnetic phases, see Fig. 3 (c). 

Turning to the spin-spin correlations, their abrupt change at 
the transition can be seen in the behavior of the ‘symmetrized’ 
spin structure factor <S(Q) shown in Fig. 3 (d), which is de¬ 
fined as 

5(Q) = ^EE e^iSgSZ) , (2) 

a r^O 

where N is the number of sites, Q i n:/ is the ordering wavevec- 
tor (see below) of the a-th component of the spins (a = 
x , y, z), and the extra factor of 2 in this definition accounts 
for the fact [33] that, for finite systems, there are no correla¬ 
tions between NN ladders like the ones shaded in Fig. 1, due to 
the non-global symmetry H xyz discussed above. These data 
show clearly the short-range (long-range) character of spin- 
spin correlations inside (outside) the QSL regions. 

This aspect can be seen more directly in Fig. 4, which 
shows the real-space spin-spin correlation profiles (Sf £“), in 
the three channels a = x,y, z, as calculated in the ground state 
of the 32-site cluster, inside the first QSL phase and slightly 
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FIG. 4. Real-space spin-spin correlation profiles evaluated at the 
ground state of the N = 32 cluster, inside the first QSL phase (ip = 
0.01-7T, left column) and inside the magnetic phase I (ip = 0.0287T, 
right column). Different rows correspond to the three different chan¬ 
nels a = x,y and z. The reference site i is indicated by the 

small black open circle. Positive (negative) correlations are shown 
by filled blue (filled red) circles, whose radius scales with the mag¬ 
nitude of the correlation. The difference between a = z and a = x, 
y stems from the fact that the 32-site cluster does not have the full 
point-group symmetry of the infinite lattice, and the momentum point 
M- is not equivalent by symmetry to M, and M a , see [33]. 


outside (magnetic phase I). The results show clearly the ultra 
short-range nature of the correlations inside the QSL region, 
and the long-range nature outside. 

Finally, the spin-spin correlation profiles demonstrate the 
special anisotropic character of the correlations, whereby dif¬ 
ferent spin components a are correlated along different di¬ 
rections of the lattice (or, equivalently, different spin compo¬ 
nents a order at different ordering wavevectors see also 
Fig. 2), reflecting the locking between spin and orbital degrees 
of freedom in this model. Similar behavior is found for all 
other magnetic phases, including the zig-zag phases that are 
relevant for Na 2 lr 03 and n -RuCL,. Such a signature of direc¬ 
tional dependent Kitaev couplings is exactly what has been 
reported recently by S. H. Chun et al. for Na 2 lrC >3 [7]; see 


also last paragraph of Sec. VII. 

In the following we shall probe the physical mechanism of 
the spin liquid instabilities by taking one step back and exam¬ 
ining the classical limit first. 


III. Classical limit 


For classical spins, the frustration introduced by the K 2 
coupling is different from the one of the pure K\ model stud¬ 
ied by Baskaran et al [47]. A straightforward classical mini¬ 
mization in momentum space [33] gives lines of energy min¬ 
ima instead of a whole branch of minima [47], suggesting a 
sub-extensive ground state manifold structure, in analogy to 
compass-like models [48] or other special frustrated antifer- 
romagnets [49], 

We can construct one class of ground states by satisfying 
one of the three types of Ising bonds. We can choose for ex¬ 
ample the horizontal zz-bonds and align the spins along the 
z-axis with relative orientations dictated by the signs of K\ 
and K 2 . The energy of the resulting configuration saturates 
the lower energy bound [33] Eb/(NS 2 ) = —I-K 2 1 — |-Ki]/2 
and is therefore one of the ground states. We can then gen¬ 
erate other ground states by noting that A'i and K> fix the 
relative signs of the spin projections S z only within the ver¬ 
tical 2-leg ladders of the lattice (shaded strips in Fig. 1), but 
do not fix the relative orientation between different ladders, 
because these couple only via xx and yy Ising interactions 
which drop out at the mean field level. This freedom leads to 
2 "“ ground states, where niadcx \fN is the number of vertical 
ladders. This sub-extensive degeneracy stems from the pres¬ 
ence of non-global, sliding operations [48, 50-52] of flipping 
S z 1 —^ - S z for all spins belonging to one vertical ladder. Sim¬ 
ilarly, we can saturate the xx or the yy bonds, leading to 2-leg 
ladders running along the diagonal directions of the lattice. In 
total, this procedure delivers 3 x 2 niad classical ground states. 

These states are actually connected in parameter space by 
valleys formed by other, continuous families of ground states 
that can be generated by global SO(3) rotations of the discrete 
states [33], The degeneracy associated with these valleys is 
accidental and can therefore be lifted by fluctuations. This is 
in fact the situation at finite T where thermal fluctuations se¬ 
lect one of the three types of discrete ground states, thereby 
breaking the three-fold symmetry of the model in the com¬ 
bined spin-orbit space. This corresponds to a finite-T nematic 
phase where spins point along one of the three cubic axes but 
still sample all of the 2 n “ corresponding states, without any 
long-range magnetic order. To achieve the latter one needs 
to spontaneously break all sliding symmetries and this cannot 
happen at finite T, according to the generalized Elitzur’s theo¬ 
rem of Batista and Nussinov [50]. The sliding symmetries can 
break spontaneously only at T = 0 and in all possible ways, 
which is reflected in the divergence of the spin stmcture factor 
along lines in momentum space. 
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IV. Quantum spins & Strong-coupling expansion 


Turning to quantum spins, the situation is fundamentally 
different because the sliding symmetries are absent from the 
beginning: To flip one component of the spin we must com¬ 
bine a 7r-rotation in spin space and the time reversal oper¬ 
ation [53]. The latter, however, involves the complex con¬ 
jugation which cannot be constrained to act locally only on 
one ladder. Essentially, this means that the ladders must cou¬ 
ple to each other dynamically by virtual quantum-mechanical 
processes, which in turn opens the possibility for long-range 
magnetic ordering even at finite T. 

The natural way to understand the dynamical coupling 
between the ladders is to perform a perturbative expan¬ 
sion around one of the three strong coupling limits where 
the above discrete states become true quantum-mechanical 
ground states. Consider for example the limit where the xx 
and yy couplings, denoted by K‘f' y> and K^ v \ are much 
smaller than the zz couplings, Kf and iff. Let us also 
parametrize = rKf 2 , Kf = cos ip and Kf = sin ip. For 

r = 0 we have W| ;l( | decoupled vertical ladders, and 2 n “ quan¬ 
tum ground states. Degenerate perturbation theory [33] then 
shows that the degeneracy is first lifted at fourth order in r via 
three, loop-four virtual processes that involve: (i) only Kf <v \ 
(ii) only Kf y \ and (iii) both Kf (v) and K% [v) perturbations, 
see the top panel of Fig. 5. 

The processes (i) give rise to intra-ladder, six-body terms 
which are nothing else than the flux operators W p . As shown 
by Kitaev [10], these terms can be mapped to the square lattice 
Toric code [54] which has a gapped spin liquid ground state. 
Next, the processes (ii) and (iii) give rise to effective, NNN 
inter-ladder couplings of the form JSfSf, where i and j have 
the same (ii) or different (iii) sublattice unit cell indices, see 
top panel of Fig. 5. To fourth-order in r, the corresponding 
couplings Jw (i), J\ (ii), and ,/ 2 (iii) read 


Jw = 


Ji = 


- (KfKfr \KI\ 


64(\Kf\+2\Ki\y(\Kf\+3\Kf\)(\Kf\+4\Kf\)’ 


(KfK. 


8(\Kf\+2\Kf\)H2\Kf\+3\Kf\) 


sgn (Kf 


J 2 = 


K f K f Kf Kf 

r \Ki\+\i<f\ 

4(|A'f | + 2|A||) 3 

[2\Kf\+3\Kf\ 


- + 


2|A'f 


Afl+4|A'f 


(3) 


Note that J 2 is always AFM and competes with J\ in the re¬ 
gions I and III of Fig. 2. We also emphasize that there is no 
SfSf coupling when i and j belong to NN ladders. This is ac¬ 
tually true to all orders in perturbation theory, because of the 
above non-global symmetry H xyz , which changes the sign of 
S, on every second vertical ladder (B and C sites of Fig. 1). 

The main panel of Fig. 5 shows the behavior of | Jw\/r A , 
2 1 Ji|/r 4 , and J 2 /r 4 as a function of the angle ip, where the 
relative factor of 2 between | .J\ | and J 2 accounts for their rel¬ 
ative contribution to the total classical energy. Close to the 
exactly solvable points ip = 0 and n, the physics is dominated 
by the flux terms W p which, as mentioned above, lead to the 
gapped Toric code QSL [10, 54], The gapless QSL at r = 1 
is eventually stabilized by off-diagonal processes that neces¬ 
sarily admix states outside the lowest manifold of the r = 0 



FIG. 5. (Color online) Top: The three types of virtual processes 
around the strong coupling limit r = 0 [33]. Bottom: \Jw\/r i , 
2|Ji|/r 4 , and J 2 /r 4 vs ip. The shaded strips denote the regions 
where J 2 competes with Ji and J 2 > 2| Ji |. 


point [55]. 

The four magnetic phases I-IV of Fig. 2 are all stabilized 
by Ji which, according to Fig. 5, is the dominant coupling 
in a wide region away from ip = 0 and n. Note that there 
are also two windows (shaded in Fig. 5) in the beginning of 
regions I and III where the two inter-ladder terms compete and 
21 Ji| < J 2 . This opens the possibility for two more states (the 
ones favored by J 2 ) in these regions. This scenario is however 
not confirmed by our ED spectra and spin structure factors 
(especially for the 32-site cluster which is commensurate with 
both types of competing phases), showing that these phases 
are eventually preempted by the QSLs and the phases I and III 
at higher values of r. 

We remark here that the 1-loop formulation of PFFRG de¬ 
livers the J 2 but not the J\ processes because, in a diagram¬ 
matic formulation of Abrikosov fermions, these processes re¬ 
late to 3-particle vertex contributions, which require a 2-loop 
formulation. However, for ip around 0 and 7r, where Ji is 
small, a 1-loop formulation already yields good agreement. 


V. Semiclassical picture 


The magnetic phases of the model can be captured by a 
standard semiclassical expansion, but this has to go beyond 
the non-interacting spin-wave level. Indeed, the zero-point 
energy of the quadratic theory lifts the accidental continuous 
degeneracy of the problem (selecting the cubic axes for the 
global direction in spin space, see Ref. [33]), but fails to lift 
the discrete 2” lad degeneracy (the spectrum has lines of zero 
modes corresponding to the soft classical twists along individ¬ 
ual ladders), and does not deliver a finite spin length, in anal¬ 
ogy to several frustrated models [31, 49, 56, 57]. The spurious 
zero modes are gapped out by spin-wave interactions, leading 
to the expected anisotropy gap and a finite spin length. The 
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latter (obtained here from a self-consistent treatment of the 
quartic theory; details will be given elsewhere) tracks closely 
the behavior of the spin length extracted from the ED ‘sym¬ 
metrized’ spin structure factor [58] <S(Q), see Fig. 3 (d). Fur¬ 
thermore, both methods give values that are very close to the 
classical value of 1/2 inside the magnetic regions, showing 
that these phases are very robust. The quartic spin wave ex¬ 
pansion is however insensitive to the proximity of the QSFs, 
most likely due to the first-order character of the transitions. 

VI. Triangular Kitaev points 

At ip = :t / the system decomposes into two inter¬ 
penetrating triangular sublattices, where the K 2 coupling 
plays the role of a NN Kitaev coupling. This problem has been 
studied for both classical [28, 29] and quantum spins [30-32], 
The above analysis for the magnetic phases still holds here, 
the only difference being that the two legs of each ladder de¬ 
couple, since they belong to different triangular sublattices. 
The ordering between the legs belonging to the same sublat¬ 
tice stems from the effective coupling J 1; which is the only 
one surviving at Ki = 0. This coupling connects NNN legs 
only, leading to twelve states in each sublattice and thus 12 2 
states in total, instead of 12 for finite K\. The accumula¬ 
tion of such extra states at low energies can be clearly seen 
in Fig. 3(a-b) at ip = ±f. Note that while the ED spectra 
are broadly independent of system size, significant differences 
between the two cluster sizes are apparent near ip = ±7r/2. 
These differences, e.g. on the ground state multiplicity, can 
be easily traced back to the different point group symmetry of 
the two clusters, see detailed explanation in [33]. 

Finally we would like to point out that the origin of the 
ordering mechanism at the triangular Kitaev points has also 
been discussed independently in a recent paper by G. Jackeli 
and A. Avella [31]. 


VII. Discussion 

Charting out the stability region of the Kitaev spin liquid 
is an extremely relevant endeavor for the synthesis and char¬ 
acterization of new materials. One of the counterintuitive re¬ 
sults of this study is that the frustrating (with respect to long- 
range magnetic order) NNN coupling K 2 , which has exactly 
the same anisotropic form and symmetry structure as the K\ 
term, destabilizes the Kitaev spin liquid much faster than the 
non-frustrating isotropic Heisenberg J\ coupling. This find¬ 
ing gives a very useful hint in the search of realistic materials 
that exhibit the Kitaev spin liquid physics. In A 2 lr 03 mate¬ 
rials, for example, the role of the size of the central ion (Na 
in Na 2 lr 03 , or Li in L^IrO.j) in mediating the K 2 coupling 
(see also below) is a key aspect that can be easily controlled 
by experimentalists [59, 60], 

On a more conceptual note, the physical mechanism under¬ 
pinning the magnetic long range ordering in the present model 
is a novel example of order-by-disorder. Unlike many other 
classical states, here the ordering manifests only for quantum 


spins and not for classical spins. This striking contrast be¬ 
tween classical and quantum spins is even more surprising in 
the light of the fact that all these phases have a strong classical 
character with local pseudo-spin lengths that are very close to 
the maximum classical value of 1 / 2 . 

On this issue, we should stress that there is no discrep¬ 
ancy between the very large pseudo-spin length that we re¬ 
port here and the small length of the magnetic moments ex¬ 
tracted from magnetic reflections, e.g., in Na 2 lr 03 [5], Such 
an apparent discrepancy can be explained by the value of the 
^-factor which can be significantly smaller then 2 , because 
the orbital angular momentum is not quenched in strong SOC 
compounds. For the ideal cubic symmetry, for example, the 
well-known Lande formula gives g = 2/3, and similar values 
could be expected for lower symmetry. 

Let us now elucidate further our main reasons on why the 
K 2 coupling must play an important role in Na 2 lrC> 3 , and can 
be relevant in L^IrCU and 0-R11CI3: 

i) The super-exchange expansion of [20] shows clearly 
that the NNN Kitaev coupling is the second largest term in 
Na 2 lrC> 3 , with K 2 ~ 7-9 meV. All other perturbations are 
at most 1-2 meV, consistent with the numbers given by the 
large-scale ab initio quantum chemistry study of [15]. The 
mechanism behind the large magnitude of K 2 in Na 2 lrOs is 
physically very clear: It originates from the large diffusive Na 
ions that reside in the middle of the exchange pathways, and 
the constructive interference of a large number of four path¬ 
ways [ 20 ]. 

In L^IrOs, the K 2 interaction comes from the same mech¬ 
anism but it is relatively smaller because of the smaller size 
of Li ions [26]. Still, as discussed in [22], this coupling can 
be important to explain the current experimental evidence in 
terms of magnetic susceptibility profile, Curie-Weiss temper¬ 
ature, and the relevant range of couplings. 

Finally, in a-RuCI>. the analogous super-exchange path is 
absent, but an appreciable K 2 still arises from the anisotropy 
of diagonal interactions originated from the interplay between 
different hopping processes [26]. However, as we already 
pointed out in the Introduction, the second largest coupling 
in a-RuCl 3 is the anisotropic exchange T [15, 27]. According 
to the study of J. Rau et al. [27], a positive T seems to compete 
with K 2 for positive K\ [26]. However, the situation is still 
unclear since the Bragg peaks of the states favored by T do 
not reside at the M points of the BZ found experimentally by 
J. A. Sears et al. [13], whereas such Bragg peaks are naturally 
present in the zig-zag phases favored by K 2 , or even by a neg¬ 
ative Ji. So a lot more work is needed to clarify the relative 
importance of T, K 2 . and J\ in Q-RUCI 3 . 

ii) The K 2 coupling explains naturally the zig-zag order¬ 
ing in Na 2 lrC> 3 . This phase cannot arise in the original J\- 
K\ model, because this would require an AFM coupling K \, 
whereas it is widely accepted that K\ is FM and large in mag¬ 
nitude, see e.g. [16]. Also, the much smaller T terms, which 
are positive, also favor the zig-zag phase and do not compete 
with K 2 , according to [27]. 

iii) The K 2 coupling can provide in addition the basis to 
resolve the long-standing puzzle of the large AFM Curie- 
Weiss temperature [2, 3, 6 ], without incorporating unrealis- 
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tically large values of longer-range Heisenberg couplings J 2 
and J 3 . 

iii) The recent diffusive x-ray scattering experiments by 
S. H. Chun et al. [7] have provided direct evidence for the 
predominant role of anisotropic, bond directional interactions 
in Na 2 Ir 03 . In conjunction with the above discussion and the 
results of Fig. 4, the K 2 term then emerges naturally as the 
number one anisotropic candidate term that can drive the zig¬ 
zag ordering and the directional dependence of the scattering 
found in [7], 

An aspect that remains to be discussed in the context of 
Na 2 lr 03 is the direction of the magnetic moments which, ac¬ 
cording to the x-ray scattering data of S. H. Chun et al. [7], do 
not point along the cubic axes but along the face diagonals. 
As discussed above, the K 2 coupling stabilizes the zig-zag 
phase but it is unable to lock the direction of the moments at 
the mean-field level due to an infinite accidental degeneracy. 
The fact that the locking along the cubic axes in the Ki-K 2 
model eventually proceeds via a quantum order-by-disorder 
process (see Ref. [33]) renders this result very susceptible to 
much smaller anisotropic interactions that can pin the direc¬ 
tion of the moments already at the mean field level. A very 
small positive anisotropic T term can for example play such a 


role and can account for the locking along the face diagonals, 
as can be directly seen by a straightforward minimization of 
the classical energy. An alternative scenario involves a com¬ 
peting order-by-disorder effect within a more extended model 
that includes weak longer-range exchange interactions [26]. 
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Supplemental material 


i 


In this Supplementing material we provide auxiliary information and technical details and derivations. Specifically, Sec. A 
deals with the Luttinger-Tisza minimization of the classical energy in momentum space (1), and the order-by-disorder process by 
harmonic spin-waves (2). Sec. B gives details about our finite-size ED study, including the symmetry analysis of the low-energy 
spectra in regions I and II of the phase diagram (3), and the definition of the ‘symmetrized’ spin structure factor <S(Q). In Sec. C 
we provide results from the pseudofermion functional renormalization group (PFFRG) approach. Finally, in Sec. D we provide 
the derivation of the effective Hamiltonian around the strong coupling limit of =0. 


A. Semiclassical analysis 
1 . Lutinger-Tisza minimization 


We choose the primitive vectors of the honeycomb lattice as ti =ay and t 2 = (—^x+^y)a, where a is a lattice constant, 

see Fig. 1 of the main paper. We also define t 3 = ti —1 2 = 2 ''x + 4y. In the following, we label the Bravais lattice vectors as 
R = nti+mt 2 , where n and m are integers. We also denote the two sites in the unit cell by a sublattice index i = 1-2. The total 
classical energy of the K \ - K 2 model reads 


E — ^2 Ei (<Sr,,i Sr .,2 + Sj 


R,1 °R+t2,2 


+ S£ 


QV 

R,1 J R-t, 


2 ) + K 2 ^2 (Sr„ i S R _ti,i + S R i Sj 


R+t3,i 


^R,i *^R+t2 ,i) 


(Al) 


Defining S R)i = ^ k e' k ' R S kti , we get 

e = H/N uc = K 1 £ [s kjl SI k , 2 + e~ ik t >SZ }1 S* k , 2 + S\ >2 ] 

k 

+ K 2 ^2 [cos(k ■ ti) Si ti S*_ Ki + cos(k ■ t 3 ) S k ,i S* kii + cos(k ■ t 2 ) SJF S“ kji ] 

k,i 

= EE^, i -A(; ) (k).s“ kJ , 

k ,ij a. 

where N UC = N/ 2, is the number of unit cells, and the matrices (where a — x, y , z) are given by 


I< 2 cos(k ■ t 3 ) At e - ik t2 


A w (k)=|^jr 2 


K 2 cos(k ■ t 3 ) 


K 2 cos(k ■ t 2 ) Ar e lk t3 


, A (a) (k)=( 2 


ikt 3 


K 2 cos(k ■ t 2 ) 


, (k)=( 

2 


l\ 2 cos(k • ti) 

K 2 cos(k ■ ti ))' 


To find the classical minimum we need to minimize the energy under the strong constraints S R j = S 2 , V(R, i). The Futtinger- 
Tisza method [1-4] amounts to relax the strong constraints with the weaker one i S R i = NS 2 , or equivalently ]T k . S k , • 



FIG. 1. The first two Brilouin zones of the honeycomb lattice, along with the special lines in momentum space Q*-^, Q iv \ and Q*- 2 ^ 
(respectively , and Q (2 ^ ) corresponding to the minima of the classical energy for K 2 > 0 (< 0), see text. 
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S_k,i = S' 2 . If we can find a minimum under the weak constraint that also satisfies the strong constraints then we have solved 
the problem. To this end, we minimize the function 

T = e-A^(S k , i -S_ k , i -S 2 ), (A2) 

k ,i 


with respect to {,S'l' k ( }, which gives a set of three eigenvalue problems for the A matrices: 


y A^f—q) Sqj = A a = x,y,z. (A3) 

3=1.2 


If we can satisfy these three relations (plus the strong constraint) with a single eigenvalue A, then e = AS 2 . So the energy 
minimum corresponds to the minimum over the three eigenvalues A <n> of the matrices A' :V A(- k), and over the whole Brillouin 
zone (BZ). The eigenvalues of these matrices and the corresponding eigenvectros are: 


a£° = AT cos(k ■ t 3 ) ± ^ AT, A^° = AT cos(k ■ t 2 ) ± * Ai, a£ z) = AT cos(k ■ ti) ± ^ AT, 


r ( x ) 


±e" 


.(») 


±e 




1 

±1 


For I \2 positive, the minima of A± \ A±' 1 , and A^T* are located on the lines = r(G-|+G 2 )+(/+;j; )G 2 , = rGi+(7-|-^)G 2 , 

and Q( z ) = r G 2 + (1+ |)Gi, respectively, where l is any integer and r £ (— |). On the other hand, for A 2 negative, the 

minima are located on the lines: Q w '=r(Gi+G 2 )+iG 2 , Q( y )'=rGi+(G 2 , and Q^' = rG 2 +ZGi. Both sets of lines are 
shown in Fig. 1. 

Let us now try to build a ground state from the minima of the above eigenvectors for the case AT, 2 > 0, by using the line of 
minima as follows: 


— S y /q(^) e 
{Q (z B 


i Q 


(*).■ 


l 

-l 


= (-1)"S 


£ m 

-G 


(A4) 


i /q 

where we used the relation R = nt 1 + mt 2 and have defined = j_ f /2 dr f (r)e l2 ’ Krnr , which is the Fourier transform of the 
envelope function fir). We still need to satisfy the spin length constraint, which imposes a condition that the inverse Fourier 
transform of f(r) takes only the values ±1. This freedom corresponds to the sliding symmetries of flipping individual vertical 
ladders, and leads to 2 n “ degenerate states (where n,| a< i is the number of vertical ladders), as discussed in the main text. 

Similarly we can construct another 2 x 2 raild states by using the lines Q lfx i or in momentum space, which correspond to 
decoupled ladders running along the diagonal directions of the lattice. Altogether, we have found the 3 x 2 raiad discrete classical 
ground states discussed in the main text by using the Luttinger-Tisza minimization method. 

Finally, it is easy to see that we can also combine the three types of states into a continuous family of other ground states that 
include coplanar and non-coplanar states. This family can be parametrized by two angles 6 and o as follows. 


S 


R ,i — 



(A5) 


where ( = 1,2 and SZi, ,Sp t and .S’fj , denote the three type of discrete solutions discussed above. 


2 . Harmonic order-by-disorder 


As we claimed in the main text, harmonic spin waves lift the accidental continuous degeneracy of the classical ground state 
manifold and select the discrete 3 x 2™ lad states, whereby spins point along the cubic axes. Here we shall demonstrate this result 
by considering a one-parameter family of coplanar states obtained by linearly combining two zigzag states and two stripy states 
with spins pointing along the cubic axes. In the resulting family of states, spins are pointing in some direction on the zx-plane. 

Figure 2 shows the two zigzag and two stripy phases with spins pointing along the cubic axes. Here “yz-zigzag//x” denotes a 
zigzag state with FM zig-zag lines running along the yy and zz bonds of the Kitaev Hamiltonian, and the spins point along the 
x-axis. Similarly, “x-stripy//z” denotes a stripy state with FM ladders formed by the xx bonds of the Kitaev Hamiltonian, and 


iii 


(a) yz-zigzag // z 


(b) zx-zigzag // x 


(c) x-stripy // z 


(d) y-stripy // x 





se/n 



5E (2) /iVuc 




8/n 


FIG. 2. Two representative zigzag (a-b) and two stripy (c-d) phases that belong to the classical ground state manifold inside the regions II and 
I, respectively. Blue (red) circles denote spins pointing up (down) along the z-axis for (a,c) or along the x-axis for (b,d). The shaded stripes 
denote the FM zigzag lines in each of the two zigzag phases (a-b), or the FM ladders in each of the two stripy phases (c-d). (e-f) Harmonic 
zero-point energy SE^ (divided by the number of unit cells N uc ) as a function of the parameter 6, for two representative points inside region 
II (e) and I (f). 


the spins point along the z-axis. Specifically, these states can be written as: 

( Sr,i \ iM x R / z \ / ,\n+m ( z \ . • „ 

\ Sr,2 J = 6 (,zj =( - 1) (z) //Z 

= e iMs ' R x ) = ( x ) “*■ zx ' zi S za S //x 

= e iMx ' R ( _ Z z ) = (-ir +m ( _ Z Z ) ^ axstripy //z 

= e iM “' R ( _ X x ) = (-1)" ( _ X x ) -*■ y-stripy//x 

where M a . = n'j and M y = f^- Tr j ( see Fig- 1) an d FI = n ^i + Tnt 2 . The one-parameter family of classical ground 

states are obtained by linear combinations of the above states: 


f Sr.i \ , / (-l) m cos6»z + sin6(x \ 

\ Sr. 2 J y C(“l) m cos ^ z + C s in^ x /’ 


(A6) 


where £ = 1 for the zigzag case and £ = — 1 for the stripy case. The effect of harmonic spin waves can be found by a 
standard linear spin-wave expansion around the corresponding states for each value of 6. Figs. 2 (e-f) show the zero-point 
energy correction (per number of unit cells) as a function of the angle 9 for a representative point inside region II {ip = 0.87T, 
£ = 1) and another point inside region l(ip = 0.37T, £ = — 1). The data show clearly that harmonic fluctuations select the states 
with the spins pointing along the cubic axes (9 = 0, ±7t/ 2, and n). 

We have checked that the result is the same for the corresponding order-by-disorder process for the one-parameter family of 
states obtained by combining two states with the same wavevector, such as the “zx-zigzag // z” and “zx-zigzag // x”. 


B. Technical details about the ED study 
1 . The symmetry group of the Hamiltonian 


The full symmetry group of the K \ -K- 2 model, for half-integer spins, is T x Cq v x D 2 , which consists of: 
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TABLE I. Character table of the double covers and D 2 of the point groups C6v and C6v, respectively. The thick horizontal line separates 
the regular from the spinor IRs. 


f - 6v 
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E 5/2 
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1 . The translation group T generated by the primitive translation vectors ti and t 2 , see Fig. 1 of the main text. 

2. The double cover C ()V of the group Cg v C S0(3) in the combined spin and real space, where the six-fold axis goes through 
one of hexagon centers. This group is generated by two operations: the six-fold rotation C (i around [111], whose spin part 
maps the components ( x , y , z) 1 —> ( y , z, x ), and the reflection plane (110) that passes through the zz-bonds of the model, 
whose spin part maps ( x, y, z) i-A (— y, — x , —z). 

3. The double cover D 2 of the point group D 2 C SO(3), which consists of three 7r-rotations C 2x , C 2y , and C 2 , in spin space. 
The first maps the spin components (x, y, z ) 1 —)• ( 2 , — y, — z), etc. 


2. Finite clusters 

In our ED study we considered two clusters with periodic boundary conditions, one with 24 and another with 32 sites, with 
spanning vectors (2ti — 4t 2 ,4ti — 2t 2 ) and (2ti — 4t 2 ,4ti), respectively. These clusters are shown in Fig. 3 (a, c). The 24-site 
cluster has the full point group symmetry of the infinite lattice, i.e. Cc lV x D 2 , whereas the 32-site cluster has the lower symmetry 
C 2v x D 2 , where C 2v contains the reflection planes (110) and (110). Turning to translational symmetry, the allowed momenta 
for each cluster are shown in Fig. 3(b, d). Both clusters accommodate the three M points of the Brillouin zone (BZ) and are 
therefore commensurate with all magnetic states of the phase diagram. The difference between the two clusters is that the three 
M points are degenerate for N = 24 but not for N = 32. 

In our ED study we have exploited: i) translations, ii) the C 2 subgroup of full Cg v point group (which is equivalent to the 
inversion I in real space through the hexagon centers), and iii) the global spin inversion which maps the local S z basis states 
It) H► It)- This operation is described by ]~[ ( erf, which is nothing else than the global 7r-rotation C 2x in spin space, divided by 



FIG. 3. The two finite clusters (with periodic boundary conditions) used in our ED study, with N = 24 (a) and 32 (c) sites, along with the 
allowed momenta in the first Brillouin zone, (b) and (d), respectively. 
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TABLE II. Transformations of the twelve states of region I under symmetry operations of the model. The phases that appear for some 
operations follow from the action of these operations on the single spin-1/2 wavefunctions |n) corresponding to the spin pointing along n or 

c“~{?x>, |y>, |z>, |x>, |y), |z}} ^ {-i|y>, e~^\z), e^x), |y>, e~ ix ^\z), -e i7r / 4 |x>}, 

■ {|x>, |y>, |z), |x), |y), |z)} {-i|x), |y), -i|z), i|x), —|y>, -i|z)}, 

(ll0) ■ (|x), |y), |z), |x), |y), |z}} n- (e-^ly), e^ 4 |x), -e^ 4 |z>, - e -^ 4 |y), -e^ 4 |x>, e -^/ 4 |z)}. 

For N = 24 and 32, the product of all these phase factors give +1. 



7ti (real sp.) 7t 2 (real sp.) C 2 (real sp.) C 2x (spin sp.) 

c 6 

(110) 

|str, ® z ) 

|str, x z ) 

|str, as z ) 

|str, x z ) 

(— l) JV/a |str, x~ z ) 

(-ir/ 2 |str,y x > 

(-l^lstr.j,-) 

|str,j/ x ) 

lstr ,y~ x ) 

str, y x ) 

1 str, y x ) 

1 str, y x ) 

(—i)-^/ 3 |str, z y ) 

(—l) JV / 2 |str, a: y ) 

|str, z y ) 

|str, z y ) 

j str, z~ y ) 

| str, z y ) 

( — l) IY / 2 |str, z~ y ) 

(~ l) JV / 2 |str, x z ) 

(-i)*/ 2 |str,0 

|str, x~ z ) 

| str, x z ) 

| str, x z ) 

| str, a:” z ) 

(—1)^*1 str, x z ) 

(-l)^ 2 |str ,y~ x ) 

( — l) iV / 2 |str, y z ) 

1 str, y~ x ) 

1 str, y x > 

1 str, t/ _x ) 

1 str, y~ x ) 

1 str, y~ x ) 

( —i) 2 ^/ 3 |str, z~ y ) 

(— 1 ) JV / 2 1 str, x y ) 

|str, z y ) 

|str, z y ) 

| str, z y ) 

|str, z y ) 

(-l) JV / 2 |sn-, 2 y ) 

( — l) JV / 2 |str, x~ z ) 

(— j) JV / 2 |str, z x ) 

|str, x y ) 

|str, x~ y ) 

| str, x~ y ) 

| str, x y ) 

(— 1) JV ’/ 2 1 str, *- y ) 

( — l) JV / 2 |str, y z ) 

( z) iV T 2 |str, y~ x ) 

str,y z ) 

1 str, y~ z ) 

str, y z ) 

1 str, y z ) 

( — l) iV / 2 |str, y~ z > 

(— 1 )^T 2 1 str, 2 X > 

( —l) JV / 2 | s tr, a: _z ) 

str, z x ) 

|str, 2 X > 

1 str, 2 _x > 

1 str, z x ) 

1 str, 2 X ) 

( — ■i) JV T 2 |str, x y ) 

(-l) A, / 2 |str, 2 - y ) 

|str, x~ y ) 

|str, x y ) 

| str, x y ) 

|str, x~ y ) 

(—l) JV ( 2 |str, x y ) 

|str ,y" z > 

( — *) JV / 2 |str, y x ) 

str,t/ _z ) 

1 str, y z ) 

1 str, ?/ _z > 

1 str, j/~ z ) 

( — l) JV / 2 |str, y z ) 

(-1 ) JV / 2 |str, 2 - x ) 

( — l) JV / 2 |str, a: z ) 

1 str, z x > 

|str,^- x > 

|str, 2 X > 

1 str, 2 x ) 

|str, 2 " x ) 

( — i) JV T 2 | s tr, x~ y ) 

(~l) IV / 2 | S tr, 2 y > 


a phase factor i N . Consequently, the energy eigenstates are labeled by: i) the momentum k, ii) the parity under C 2 (‘e’ for even, 
‘o’ for odd), and iii) the parity under S z spin inversion (‘Sze’ for even, ‘Szo’ for odd). 


3. Symmetry spectroscopy of classical phases 

Here we derive the symmetry decomposition of the twelve magnetic states of region I and II of the phase diagram. As 
explained in the main paper, the other two regions. III and IV, map to I and II, respectively, by the hidden duality of H yxz 
followed by a simultaneous change of sign in K\ and /\ 2 . 


a. Phase 1 

In the following, | str, a^) denotes the stripy state with FM ladders running along the direction of the o-bonds, and the spins 
pointing along (3 in spin space. The twelve magnetic states of region I of the phase diagram can be split into four groups: 

51 = {|str, a: z }, |str, y x ), |str, 2 y }}, Si = {|str,:r" z ), \str,y~ x ), |str, 2 _y }}, 

5 2 = {|str, y z ), |str, z x ), |str, x y }}, S 2 = {|str, y~ z ), |str, z~ x ), |str, aT y }}, 

Table II shows how these twelve states transform under some of the symmetry operations of the group. Let us first examine the 
translation group. We have, V/3: 

7ti • | str, a;* 3 ) = | str, x~ 0 ), 7t 2 • | str, x?) = |str, x _/3 ), 

Til ■ l str , y 13 ) = |str, y~^), 7t 2 • lstr,?/ 3 } = lstr,/*), 

7ti • | str, z?) = | str, z 13 ), 7t 2 • | str, z 13 ) = |str, z~ P ). 

Thus ^ (|str, a;^) + |str, x~^}) transforms as k = 0 (T point) and (|str, x&) — | str, x~^}) transforms as k= T(— ~^,n) = 
Mj. Similarly, (| str, j/^) + |str, y~^)) transforms ask = 0and J- (|str, y^) — |str, y^ 13 }) transforms as k= 7r) = M y , 

^ (|str, z /3 ) + |str, z~P)) transforms as k = 0, and (|str, z&) — |str, z~P)) transforms as k=T(^,0) = M-. Altogether: 

{| str, x z ), | str, aT z }} -> V ® M„ {|str, y x ), |str, y _x )} -> T ® M„, {| str, 2 y ), |str, z~ y )} -> F ® M„ 

{|str, x~ y ), |str, * y )} -¥ T ® M„ {|str, y~ z ), |str, y _z }} -A T ® M„ {|str, z _x ), |str, 2 X }} -> T © M,. 

Next, let us examine the parities with respect to the C 2 rotation in real space and the C 2x rotation in spin space. It is easy to 
see that the first symmetry is not broken by any of the twelve states, while the second is broken when (3 = y and z. So all twelve 
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states are even with respect to C 2 , the (3 = x are even with respect to C->x, while (3 = y and z must decompose into both even 
and odd parities with respect to C 2a; ■ Altogether: 

{|str, x z ), |str, £ -z }} —►r.e.Sze © M^.e.Szo, {|str, x~ y ), |str, x y )} —>T.e.Sze © M :r .e.Szo, 

{|str, y x ), |str, y _x }}—►r.e.Szeffi M y .e.Sze, {|str, |str, y 2 ')} —>T.e.Sze © M y .e.Szo, (Bl) 

{|str, z y ), |str, z~ y )} —►T.e.Sze © M-.e.Szo, {|str, z _x ), |str, z x }} —>T.e.Sze © IVU.e.Sze . 

‘Extra’ degeneracy at the M points for N = 24. The above quantum numbers for the M points are fully consistent with what 
we find in the low-energy spectra of Fig. 3 (a) of the main paper. For the symmetric, TV = 24 cluster, the three M points are 
degenerate due to the six-fold symmetry. However we see that the two sets of M points are also degenerate with respect to each 
other, i.e. we have a six-fold degeneracy. This extra degeneracy comes from the D 2 symmetry in spin space. To see this, let us 
relabel the spin inversion part of (Bl) using the actual IR of the group D 2 (see Table I, right), instead of the parity with respect 
to C 2 „. (which contains less information about the state): 

{|str, x z ), |str, z )} — I 17.e.A © M^.e.Bi, {|str, x~ y ), |str, x y )} — >T.e.A © M^.e.Ba, 

{|str, y x ), |str, t/ _x }}->-T.e.A © M y .e.B 3 , {|str, y~ z ), |str, y z )} —^T.e.A © M y .e.Bi, (B2) 

{|str, z y ), |str, z~ y )} —► I’.e.A © M z .e.B 2 , {|str, z~ x ), |str, z x )}—>r.e.A © M-.e.B 3 . 

We see that the two states belonging to a given M point transform differently under D 2 , so the Hamiltonian does not couple the 
two states. Yet, these states are mapped to each other by one of the reflection planes of C6 V , so they must be degenerate, leading 
to an overall six-fold degeneracy at the M points. 

Degeneracies at the T point for N = 24. The little group of the F point is the full point group Cf, v x D 2 . However, all of the 
above six states that belong to the T point belong to the identity IR of D 2 , so it is enough to decompose them with respect to the 
Cqv part of the little group. To this end we use the well known formula from group theory [5] 

m a = -J— X a {g)X{g)* , (B3) 

|Qv| se c 6v 

which gives the number of times m a that the a-th IR of Cg v appears in the decomposition of the 6x6 representation formed by 
the six states belonging to the T point. Here A'(g) gives the character of this representation, while x“(g) is the character of the 
a-th IR of C6v, see Table I (left). From Table II it follows that A'(g) is finite only for the elements E, E, C 2 , and C 2 , and using 
the characters of Table I (left) we find that the only finite m a are the following: = toa 2 = 1, wis 2 = 2, namely 

6r -> Ai © A 2 © 2E 2 . (B4) 

i.e. we expect two singlets and two doublets. All states are found in the low-energy spectra shown in Fig. 3 (a) of the main paper, 
where the degeneracy of the E 2 levels has been confirmed numerically. 


b. Phase II 

Here we denote by zig. ao 1 ^) the zigzag state with FM lines formed by consecutive a and a’ type of bonds, and the spins 
pointing along (3 in spin space. The twelve magnetic states of region II can be split into four groups: 

5 3 = (|zig,yz z ), |zig ,zx x ), |zig, xy y )}, S 3 = {|zig, yz~*), |zig,zaT x ), |zig, xy~ y )}, 

5 4 = {|zig, zx z ), |zig, xy x ), |zig, yz y )}, S 4 = {|zig, zx _z ), |zig, xy~ x ), |zig, yz~ y )}, 

Under T and C 2 .„ in spin space, these states transform in analogous way with the twelve states of region I, see (Bl). The 
difference is that the present states break the C 2 rotation around the hexagon centers, and therefore the decomposition will 
contain both even and odd parities with respect to C 2 . Specifically, 

{|zig, yz z ), |zig, 2 /z -z }} —^T.e.Sze © M^.o.Szo, {|zig ,yz~ y ), |zig, yz y )}^> T.e.Sze © M^.o.Szo, 

{|zig, zx x ), |zig, z* _x )} —>-r.e.Sze © M y .o.Sze, {|zig, zx~ z ), |zig, z* z )} —>T.e.Sze © M y .o.Szo, (B5) 

{|zig, xy y ), |zig, xy~ y )}— >■ T.e.Sze © M-.o.Szo, {|zig, xy~ x ), |zig, x y x }}—>■ T.e.Sze © M 2 .o.Sze . 

In analogy with region I, for the symmetric 24-site cluster, the six states belonging to the M points are degenerate due to the 
additional D 2 symmetry, and the six states belonging to the T point decompose as in (B4), namely 6T —► A x © A 2 © 2E 2 . Again, 
all states are found in the low-energy spectra shown in Fig. 3 (a) of the main paper. 



c. Special points 0 = ±7r/2: Different ground state structure for N = 24 and N = 32 
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As shown in Figs. 3(a) and (b) of the main text, the ED results are broadly independent of system size but significant differences 
between the two cluster sizes are apparent for the GS structure near 0 = ± 7 r/ 2 . The reason behind this difference lies in the 
different point group symmetry of the two clusters. The 24-site cluster has the full point group symmetry of the infinite lattice, 
whereas the 32-site cluster does not. This is also true for the two triangular sublattices of each cluster at 0 = ± 7 r/ 2 , where they 
become independent from each other. 

Due to the high symmetry, each of the 12-site sublattices of the 24-site cluster have a two-fold degenerate ground state at 
0 = ±7t/ 2; let us denote them by |a) and |/3). On the other hand, the lower symmetry of the 16-site sublattices of the 32-site 
cluster leads to a single, non-degenerate ground state; let us denote it by I 7 ). Now, the global ground state structure of the two 
clusters at 0 = ±7 t/ 2 follows simply by taking the tensor product of the ground state manifolds in each sublattice. The 24-site 
cluster has four ground states: 

|ct)subl ® |ct)sub2; |/3)subl ® |/^)sub2? |t0subl ® |/^)sub2i |/3)subl ® |ct)sub2- (B 6 ) 

The first two states belong to the representation T.e.Sze, i.e. they have even parity with respect to inversion through the middle 
of the hexagons (this operation maps one sublattice to the other), and the same is true for the combination ( |a) su bi <8> |/3) su b2 + 

|/3)subi ® |c0 S ub2)- The remaining, antisymmetric combination, 4j(|a) su bi ® |/S)sub 2 - |/3)subt ® |«)sub 2 ), belongs to T.o.Sze, i.e. 
it has odd parity. This is in perfect agreement with the ED data. 

For the 32-site cluster on the other hand, there is only one global ground state, namely |'y)subi <S> | 7 )sub 2 , which has even parity, 
again in agreement with the ED data. 

Of course, as we discuss in the main text, in the thermodynamic limit a large number of states (12 2 ) will collapse to the ground 
state, which is how the corresponding symmetry-broken (classical) states are eventually formed. 


4. ‘Symmetrized’ spin structure factor and spin length 


Here we discuss the ‘symmetrized’ spin structure factor S( Q) and explain the overall normalization factor that we use to 
extract the spin length. As we discuss in the main text, NN ladders do not couple by the symmetry H xyz , and so the quantum 
ground state of a finite cluster contains both relative orientations of the two sets of ladders Li and L > with equal amplitude. 
As a result, the spin-spin correlations between two spins that belong to L\ and L 2 are zero for any finite cluster. If we wish to 
calculate the local spin lengths from the ground state spin-spin correlation data we can calculate the ‘symmetrized’ spin structure 
factor for one of the two subsets of ladders only, say L \: 

= E (5“5“,)e <Q( “ ) - ( '- r ' ) > (B7) 

1 ot r,r'Gl/i 

where iVi = N/2 is the number of sites inside the sublattice L \, and Q i a) is the ordering wavevector corresponding to the spin 
component a = {x, y. z ) . By translation symmetry, 

(S?S?) = (S? +s S?, +s ) => 5r(Q) = }EE <W>e <Q(a E (B8) 

1 a reii 

where we have chosen a reference site r' = 0. The local spin length m is then given by m 2 = -^<Si(Q). 

By contrast, the corresponding ‘symmetrized’ spin structure factor of the full lattice <S(Q), defined by 

5 (Q)=i^E E <^S r > lQ(Q) - (l - r '\ (B9) 

a r,r' 6 L1UL2 


would give in the present case 


S(Q) = ^t(Q), 


(BIO) 


and the corresponding local spin lengths would be off by a factor of \/2. 
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FIG. 4. Static spin-structure factor x“(k) plotted in the extended Brillouin zone (black lines inside the plotted region mark the boundaries 
of first Brillouin zone) for various values of 0 in the Kitaev spin-liquid phase. Note that x x:r (k) (x yy (k)) are related to x 2Z (k) by clockwise 
(counterclockwise) 27r/3-rotations in k-space. 


C. Pseudofermion functional renormalization group (PFFRG) approach 

In addition to ED we studied the K 1 -K 2 honeycomb model using the pseudofermion functional renormalization group (PF¬ 
FRG) approach. Rewriting the spin operators in terms of Abrikosov auxiliary fermions, the resulting fermionic model can be 
efficiently treated using a one loop functional renormalization group procedure. This technique calculates diagrammatic contri¬ 
butions to the spin-spin correlation function in infinite order in the exchange couplings, including terms in different interaction 
channels: The inclusion of direct particle-hole terms insures the correct treatment of the large spin limit S — > 00 while the 
crossed particle-hole and particle-particle terms lead to exact results in the large N limit. This allows to study the competition 
between magnetic order tendencies and quantum fluctuations in an unbiased way. For details we refer to reader to Ref. [6]. 

The PFFRG method calculates the static spin-structure factor as given by 

/>00 

X a/3 (k)= / dr(T T {5 a (-k ) r)S^(k ) 0)}), (Cl) 

Jo 

with 

S“(k, T ) = -!= Y, e~ ikri e HT S*e- HT , (C2) 

where r denotes the imaginary time and T t is the corresponding time-ordering operator. Being able to treat large system 
sizes (calculations for the K 1 -K 2 model are performed for a spin cluster with 265 sites) the PFFRG yields results close to the 
thermodynamic limit. Fig. 4 shows three representative plots for the momentum resolved spin-structure factor x 22 (k) in the 
Kitaev spin-liquid phase in the vicinity of 0 = 0. While in the exact Kitaev limit 0 = 0 the PFFRG reproduces the well-known 
nearest neighbor correlations as indicated by a single harmonics profile of the spin-structure factor, deviations from '0 = 0 lead 
to longer-range correlations and a more diverse spin-structure factor. 


D. Strong-coupling expansion 


Here we provide some technical details on the derivation of the effective model around the strong coupling limit of K f (y> = 
K^ y) = 0. In this limit we have «| a d decoupled vertical ladders (which are the ladders made of the z^-bonds), leading to a 
sub-extensive ground state degeneracy. The ordering pattern within each individual vertical ladder is fixed (up to a global sign) 
by the signs of K 2 and Kp. The GS degeneracy is lifted by the transverse perturbations K : f v> and K^ v \ which give rise 
to effective couplings between the ladders (or more accurately between NNN ladders, as discussed in the main paper). These 
couplings can be found by degenerate perturbation theory. Fet us denote by Hq the sum of all K f interactions and by V the sum 
of all remaining terms of the model. In the following we define the strong coupling parameter r to be the ratio between icfif 1 
and /if 2 - as in the main text. 

To write down the effective Hamiltonian we should define the corresponding Hilbert space on which it acts. Obviously, for 
Kp 0 0 this is the ground state manifold of H$, namely the 2™“ states corresponding to all possible relative orientations of the 
vertical ladders. However, special care must be taken at Kp = 0 where different rungs of a given vertical ladder do not interact 
with each other and the relevant Hilbert space is enlarged from 2 niad to 2 A/ / 2 , where N is the number of sites. To treat both 
Kp 00 and Kp = 0 cases at once we must then take the enlarged manifold of 2 ; ' v / 2 states. 
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FIG. 5. Effective couplings arising from only NN perturbations (A), only NNN perturbations (B), and both NN and NNN perturbations (C-H). 


With this definition of the target Hilbert space, the very first term of the effective Hamiltonian is a first-order coupling between 
the rungs which is proportional to iff, which fixes (for iff ^ 0 ) the relative orientation of different rungs within each vertical 
ladder. It is easy to see that the remaining degeneracy between different ladders is lifted in fourth-order in V. The effective 
Hamiltonian (up to fourth order) is then described by the expression 

Heff = Hi + PVRVRVRVP , (Dl) 

where “Hi contains the iff terms, P is the projection into the enlarged manifold of 2 iv / 2 states discussed above, and R= J 0 l I fj a 
is the resolvent, where E 0 = ( | K\ |/2 — \K 2 \)N is the ground state energy at V = 0. By expanding the different terms of V 
in (Dl) we get three types of loop-four virtual processes, that involve: i) only NN perturbations K'f vl (Sec. 1), ii) only NNN 
perturbations K^^' 1 (Sec. 2), and iii) both K'f v> and perturbations (Sec. 3). 


1. Effective terms arising from K* (y) only (Toric code terms) 


The if^ y) perturbations give rise to intra-ladder, six-body terms of the form J\yW p , where W p is Kitaev’s [7] flux operator: 

W p = 2 6 SrS'fSfS'|S'fS'| , (D2) 

where 1-6 label clockwise the six sites of the hexagon plaquette p, as shown in Fig. 5 (A). To find ./ M - in fourth order in r, it 
suffices to consider one hexagon only. Let us denote the local configuration of this hexagon in any of the ground states at r = 0 
by |fji, 02 , 03 , 0 - 4 , as, ae), with the spin projections S? = \<Jj, and aiag = a 3 a 4 = —sgn(/ff). In this case, the perturbation 
V = A can be written as (see Fig. 5): 


A = Aa + A b + A c + A d = KfSfSg + + tffSfST + , 


(D3) 


and Eq. (Dl) contains 24 terms in total, which have the form 


H (A,dcba) = p AdRAcRAbRAa p ; etC. 












X 


In the following we define n = (KfKf ) 2 and use the relations S x \a) = \ \ — a) and S y \a) = \ — a). The energy excitations of 

various intermediate states are 


A12 — Ai6 — A34 — A45 — — \Kl\ — 2IA2I, 

A 26 = A 35 = HA'f| - | K% |, 

Al 635 = A1235 = A2634 = A2645 = —|Aj | — 3 |Af|, 
Al 634 = A1245 = — 2|Af I — 4IA2I, 

A1234 = Al 645 = — I ATi | — 4 | A2 I, 


Let us first consider the terms of the type 


/t / ( A,dcba ) 


101, 02, 03,04, 05, 06} 


/I 01020405 

4 4 A45 A35 A1235 


I01, 02, —03,04, —05, 


-06} ■ 


The final state is not the same as the initial one, but belongs to the enlarged manifold of 2 ,v / 2 states, so this is a valid process. 
The operator that does the job is: 


4 A 45 A 35 Ai 


CV Q x Q z QV C x — 
"01 0 2 03 04 *J 5 a 6 — 


- - - Wr, . 

2 8 A 45 A35A1235 


This result can be also found right away by taking 


a/' 

/T-p, 


( A,dcba ) sjj (A,abed) 

— ft, 


_ ..JD Q x Q x JD QV QV JD Q x Q x JD qy QV JD V r* ( Q x QV \ QV Q x ( Q x QV \ QV C x _ 

— / 1 -r D 1 Oq 1 X 0 1 0 2 4 1,^3 *~ > 4 °5 r ‘* -jy V°1 D 1 ) °2 °3 l°4 *~>4 ) °5 °6 — — 




with Z? 1 = A 4 5 A 35 A 1 235 . Similarly 

j(A,cdba) sjj(A,abdc) ^/(A,dca6) sjj(A,bacd) ^ ,(A,cdab) , (A,bade) ttt 

^•eff — ™eff — "eff — ™eff ~ — "eff — — ™eff — 28 ]J 1 ^ P ’ 


where we used A 3 4 A 3 5 A 12 35 = A 4 5 A 35 A 12 35 = A 4 5 A 35 A 16 3 5 = A 45 A 35 A 1.235 = A 34 A 35 A 16 35 = A 45 A 35 A 12 35 = D l . So 
the eight processes {abed, bacd, abdc, bade } and {deba, dcab, cdba, edab} cancel each other out. 

Next come the processes: 


K 


( A,acbd ) 


=K 


(A,dbca) 


_ .. p QV QV JD QV qV jd q x q x jd q x q x jd 

— fli ±x,0^ 0 2 ilO3 1\01 oq 1 


D 2 


I ay c x \ <?y c x r qv a x \ ay a x _ e- 
I s ! *41 ) 02 03 \04 04 ) 05 06 — “ 28 D 


w„ 


-is; 


-is; 


with D 2 = Ai 6 Ai 634 A 2 634. Similarly, 

/t j(A,cabd) rtj(A,dbac) ^ j(A,acdb) ^_/(A,6dca) ^ j(A,cadb) qj(A,bdac ) rxr 

n eff ~ n eff ~ ™eff ~ H eff ~ ^eff ~ ^eff ~~ 2 8 D 2 


where we used A 34 Ai 6 34 A 263 4 = Ai 6 Ai 6 34 Ai 6 35 = A 3 4Ai 63 4Ai635 = D 2 . These eight processes {cabd, acbd, cadb , aedb} 
and {dbac, dbca, bdac, bdea} give the same contribution and, thus, do not cancel out. 

Finally, there are the processes 


n_/(A,cbad) _ sjj(A,dabc) _ JD QV QV D Q x C x JD QV QV p Q x Q x JD 

ft e ff —/L e ff —/xjT O5 iiJi Oq JADi 02 ± 0 , 0 ^ O4 jT - 


d 3 


(sr s{) s y 2 s x 3 (s y s:) s y 5 s%=+ 


2 s D 3 ^ p ’ 


is; 


-is; 


with -D 3 = A 34 Ai 234 A 2 g 34 . Similarly 

j(A,bcad) _ sjj(A,dacb )_ sjj(A,cbda) _ ^j{A,adbc) _ ^j{A,bcda) _ sjj(A,adcb) _ . ttt 

~ H eff ~ n eff ~ H eff ~ ^eff — H eff ~ + 2 8 D 3 ^ P ’ 

where we used Ai 2 Ai 23 4 A 26 34 = A 3 4 Ai 2 3 4 Ai 2 35 = Ai 2 Ai 2 34 Ai 23 5 = D 3 . So the eight processes {ebad, bead, cbda , beda} 
and {dabc, dacb, adbc, adeb} also do not cancel out. Altogether: 


nj(A) onj(A,acbd) 

n eff ~ 5 n eff 


+ 81 -L 


( A,cbad ) 
eff 


iit_L 

2 s ' D 3 


— \w — M(Al 634 — A1234) 

D 2 P 2 5 A 34 A1634 A1635 A1234 P 


We have A1634 — Ai 2 34 = — |ATf |, and therefore 


T „,W 


r 1 

'‘'eff D w vv p 


w 26(1^1 + 2\K;\)*(\Ki\ + 3|A||)(|A'f I + 4|A'||) 


(K x K y ) 2 

For A'l = 0 we get J\y = — 2 g\k‘\ 3 • w hich agrees with the result obtained by Kitaev [7], 
















XI 


2. Effective terms arising from K{ {y) only. 


Consider three consecutive ladders in the honeycomb lattice. We will show that the terms give rise to an effective NNN 

inter-ladder coupling of the form J\ S{S{, see Fig. 5 (B). In this case, the perturbation V = B is given by (see Fig. 5): 


B = B a + B b + B c + B d = K%SiS% + + Af S?S? + . (D5) 

Again, Eq. (Dl) gives 24 relevant contributions. In the following we define A = (A'fA'|) 2 , and use the relation 173(75 = 
—sgn(Af). We also introduce the excitation energies of various intermediate virtual states: 

A13 = A17 = A15 = A37 = A57 = — |Af| — 2 |Aj|, A35 = — \K{\ — |A||, Ai357 = — 2 |Af| — 3 |A|| . 


We find: 


A, 


(B ,abcd) I 


A, 


(B ,abdc) I 


Cl,(72,(73,(74,(75,(76,77} = + 
71,72,73,74,75,(76,77} = — 


Act 1(77 


4 4 Ai3Ai7Al5 


— Sgn(A'f )| 7 l, 72 , 73, 7 4 , 75, 76,77} 


Act 7(71 


4 4 Ai3Ai7A 


n, 


(B,bacd) I 


71,72,73,74,75,76,77} = S= — 


I34A17ZA57 

A7177 


4 4 A37Al7Al5 


-Sgn(A'|)|7l, 72, 73,7 4 , 75, 76,77) , 

Sgn(A'!)| 7 l, 7 2 , 7 3 , 74 , 75 , 76 , 77 ) , 


A ' badc) \cr 1 , 72 ,73, 74, 75, 76 , 77} = + . . Ag ^ <7 ' . -Sgn(Af )| 71 , 72 , 73, 74 , 75, 76 ,77} . 

4^ZA37ZAi7ZA57 

So the eight terms coming from {abed, abdc,bacd, bade] cancel out, and the same is true for their inverse processes 
{deba, cdba, dcab, edab}. Next come the processes: 

,c 6 da) |7l, 72 , 73 , 74 , 7 5 ,7 6 , 7 7 ) == ~ , 4 ,^' a -Sgn(A'|)1 7l , 72 , 7 3 , 7 4 , 7 5 , 7 6 , 7 7 ) , 

ZA13ZA35ZA57 

A e ( f f’ C6a ‘ 2) |7l,72,73,74,75,76,77) = = + a -Sgn(A'|) 17l, 7 2 , 7 3 , 7 4 , 7 5 ,7 6 , 7 7 ) , 

4 4 Ai5Z\35A57 

and similarly H^' bcad> = —H^’ cbad \ and H^’ bcda ' > = —H^’ cbda \ So the processes coming from {ebad, cbda, bead, beda} 
cancel out, and the same is true for their inverse processes {dabc, adbc , dacb, adeb}. 

The only finite contributions then come from the remaining eight processes: {acbd,acdb,bdac,bdca} and their inverses 
{dbca,bdca,cadb,acdb}. Here H^’ acbd ' > = H^’ cabd ^ =H^’ acdb ' > = H^’ cadb \ so there is no cancellation. We have: 

A(S,dtea)|7i, 7 2 ,73, ^ (j 5 , ( Jq , 77} = - XrJ l a ‘ -Sgn(A 2z )| 7 i, 7 2 , 7 3 , 74 , 7 5 , 7 6 , 7 7 ) . 

4^ZAi3Ai357ZXi5 


In total, the effective terms arising from the NNN perturbations K^ v> is 


--7 /(B) _ q / tj (B,dbca') _ j qz qz 

n eff ~ — JlDi^7 



(D6) 


For Kf = 0, Ji = 24 - 4 (a : 2 -)3 s g n (ATf), in agreement with the result obtained by Jackeli and Avella [8] for the triangular lattice 
case. 


3. Effective terms arising from mixed K{ ly> and K{ iy> perturbations. 

Finally, we consider the perturbations due to mixed and l({ ‘ y 1 terms. Figure 5 (C-H) shows the six minimal loops that 

contribute to an effective coupling of the form J 2 S'j : S'|, between sites 1 and 4. In the following we define n = Kf KfK%K %, 
and introduce the excitation energies of various intermediate virtual states: 

A12 = Aie = A14 = A46 = A24 = — | Af | — 2 |A 2 |, 

A 2 6=A35 = — |Ai| — |A 2 |, A 2 3 — A56 — 2 | A 2 [, 

Al 246 = A1345 =— 21 K{ I — 3| Aj|, Al 2 34 = A1456 = —\K{\ — 4|A||. 

Let us discuss the different processes C-H of Fig. 5 one by one. 














a. C & D processes 


The perturbation V = C described by the loops of type C of Fig. 5 splits as 

C = C a + C b + C c + C d = K}S}S% + K v 2 S v 2 Si + KZStSZ + K}S^S{. 
Replacing (D7) into (Dl), we get twenty four contributions. We have 


(D7) 




|(7l,(72,<73,<74,<75,<76) = C |(7i,<72,(73,(74,(75,(76) = “ 


4 4 Al2Al4Ai6 


I <71, 02, 03, (74, 05,(76} . 


We also find li!^ Acah> = 'H^’ cdba ' 1 ='H < ^’ dcha \ So all eight processes {dcba,dcab,cdba,cdab} and {abed,bcicd,abdc,bade} give 
the same contribution. Next come the processes of the type 


n AC.dbca) \ \ n AC.acbd) i \ 

rf e ff pi, (72, ( 73 , cr 4 , cr 5 , a 6 / — rt c ff |<7i, <72, <73, <7 4 , <75, <7 6 ) — 


K0104 


4 4 A12 A1246 A13 


|<7l, 02, CT3, O4, <7 5 , <7 6 ) 


and 'Hgg’ dbac ^ = 'H < ^’ bdca ' > = 'U^' dbca \ So all 8 processes {dbca, dbac,bdca,bdac} and {acbd, cabd,acdb,cadb} give the same 
contribution. Finally there are the processes of the type: 


0/1 


I \ n AC,adbc) < \ 

\Ol, 02, (73 , 04, (75, <76/ = H eff |<7i, (72, (73, (74, (75, (76/ = —- 


|<71, 02,03, 04, 05,06) . 


4 4 A 12 A 26 A 46 

Flere, however, -}{^ cbad '> = —%^ ,cbda \ and similarly j{d^' bcda '> = _^(C,cbda). ^ a resu |^ the last eight processes 

{cbda,bcda,cbad,bcad} and {adbc,adcb,dabc,dacb} cancel out. So the total contribution from the C loops of Fig. 5 (C) is: 


n_l(C) onJ (C,abcd) . q-u (C,dbca) K (Al 2 — Al 246 ) 
rt cff — Ort cff + 8rt ef f — oo A 3 A <71(74, 


32 A?,A 


12^ 1246 


where A 12 — A 1246 = K{ | + \K$ \ > 0. So the coupling is AFM. 
Finally, by symmetry, ='H^p. 


b. E & F processes 


These processes give rise to an overall constant, so they can be ignored. 


c. G & H processes 


Here the corresponding perturbation can be written as 

G = G a + G b + G c + Gd = K}S\S 2 + K 2 S 2 Sf + KfSi St + Af SfSS ■ 

We have 


(D8) 




(G,dcba) I 
eff 


\ n/ (G,abcd)\ \ . 

(Tl, ( 72 , ( 73 , CF 4 , ( 75 , (76/ = H ef{ \ ( 71 , (72 , (73 , O’4 , & 5 ,06/ = ^ 4 ^ 3 —|(7 1 , ( 72 , & 3 , (74, & 5 , (76/ 


—/^(7if74 . 


where we used the relation A 13 = A 14 = A 12 . Similarly, we can also show that dcba ■* =T-L^' cdbad> = 'H^’ dcab \ so the eight 
processes {dcba,dcab,cdba,cdab} and {abed,bacd,abdc,bade] give the same contribution. Next come the processes of the type: 


O/V 


I \ n AG,acbd) 1 \ 

|(7l, (72, (73, (74, (75, (76/ = rt e ff I a 1 > j (73 5 (74, CT 5 , (76/ = 


4 4 Ai 2 A1234A13 


|<7l, (72,(73,(74, (75, (7e) • 


Again, %Fp dfcca ) _^(G,dbac) _y^(G,bdca) a p e jgjj t p roce sses {dbca,dbac,bdca,bdac} and {acbd,cabd,acdb,cadb} give the 
same contribution. Finally there are the processes of the type: 


r>,(G,cbda)\ \ n/ (G,adbc)\ \ 

T^eff |(7l, (72, &3, (74, (75, 06/ = tt eff |(71 ,(72,(73,(74,(75,(76/ = — 


K,(J 1 (7 4 


4 4 A12 A23 A34 


|<7l, (72, (73, (74, <75, (76) • 


Similarly, ' cbad ' 1 ='H < ^ ,bcda ' > = — ’ cbda \ So here, {cbda,chad,beda,bead} and {adbc,dabc,adcb,dacb} cancel out. 
Altogether 

u(G) _ a Tj(G,dcba) onj (G ,dbca) _ « (Al 2 — A1234) 
n. eff — o/l cff +8/t rff — o n A 3 A <71(74, 


32 A?o A1234 


where A 12 — Ai 234 = 2|/\|| >0. So is also AFM. Finally, by symmetry, Td^p = Tdpp. 











d. 


Final result 


xiii 


nj(C—H) _ c)nj(C) . r\nj(G) _ t qz qz 

/L e ff — ^ /t e ff "r" ^ ^eff — J2*~5l^4 


J, K 

r \Kt\ + \m 2 \ki\ 1 

2 4A ? 2 

[2\Kl\ + 3\KI\ \Kf \ +4|iT||J 


(D9) 
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